Bloch function description of a Bose-Einstein condensate in a finite optical lattice 



00 
ON 



o 
o 



I. INTRODUCTION 



M. J. Steel and Weiping Zhang 
School of Mathematics, Physics, Computing and Electronics, Macquarie University, North Ryde, New South Wales 2109, 

Australia. 
(February 1, 2008) 

We consider stationary and propagating solutions for a Bose-Einstein condensate in a periodic 
optical potential with an additional confining optical or magnetic potential. Using an effective 
mass approximation we express the condensate wavefunction in terms of slowly-varying envelopes 
modulating the Bloch modes of the optical lattice. In the limit of a weak nonlinearity, we derive a 
nonlinear Schrodinger equation for propagation of the envelope function which does not contain the 
rapid oscillation of the lattice. We then consider the ground state solutions in detail in the regime 
of weak, moderate and strong nonlinear interactions. We describe the form of solution which is 
ON ' appropriate in each regime, and place careful limits on the validity of each type of solution. 
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There have now been a large number of experiments and an enormous quantity of theory devoted to the properties of 
the trapped Bose-Einstein condensates (BEC) Q|. Motivated by the existing experiments, the bulk of the theoretical 
work so far assumes a harmonic trap in one, two or three dimensions. Lately, a few papers have considered the 
I ' consequences of exposing a BEC to an optical periodic potential [|]-|| • While not yet observed experimentally, such a 
situation could easily be produced by applying counter-propagating optical beams along one or more axes to produce 
t^J- ■ a standing-wave optical lattice. Indeed, the directional output coupler of the BEC demonstrated by Leng et al. || 
uses such counter-propagating beams to induce Raman transitions into untrapped states. In this paper, we consider 
the special case of a BEC confined in a finite optical lattice, and develop theory to describe the properties of the 
O ■ lattice BEC. 

Subjecting a BEC to a periodic potential opens the way for a number of new phenomena. S0renberg and M0lmer 
assumed an infinite lattice potential and showed as one would expect, that the quasi-particle excitation spectrum 
exhibits a band-structure. Zobay et al. Q considered a propagating field of three internal levels connected by Raman 
■ transitions. They showed that if the center wavenumber is tuned close to the band gap induced by the optical lattice, 
one can derive a pair of coupled mode equations for forward and backward propagating waves. These equations 
I i are well-known to support the interesting class of solitary waves known as "gap solitons" @||. Gap solitons have 
been observed when a periodic structure in an optical fiber is subjected to high intensity light [fij but have not been 
observed in an atom optics context. Finally, Jaksch et al. Q have demonstrated an equivalence between the system 
of a condensate trapped in an optical lattice and the Bose-Hubbard model of condensed-matter physics. At very 
low site-occupation numbers, they predict the occurrence of phase transitions between a Mott insulating state and 
coherent superfluid flow. Jaksch et al. also considered a complication not present in the other two papers — the 
inclusion of a weak harmonic trap in addition to the rapidly oscillating lattice. 

It is this complication that we pursue in the present paper. As the eigenstates or "Bloch functions" of a pure lattice 
potential are of infinite extent, any real system must always include some form of additional confining potential. (In 
any case, a real optical lattice is always finite). Thus we consider the case of an infinite periodic potential combined 
with a weak confining potential. While we will allow the form of this potential to be rather general, in many cases 
a harmonic oscillator form may be appropriate. For example, one might use optical beams to generate the lattice 
and provide the ultimate confinement with a weak magnetic trap in the standard way. In fact, red-detuned counter- 
propagating focussed Gaussian beams would effectively create both a lattice and a harmonic-like axial trap due to 
the reduced intensity away from the beam focus. In this case, however, the lattice potential would also be modulated 
by a Gaussian envelope. We shall restrict to the case where the periodic part of the potential is uniform. 

Our aim is to describe in detail the properties of the condensate wave function in such a potential in both stationary 
and propagating regimes. We assume the condensate may be described by the standard Gross-Pitaevskii equation 
(GPE) with the addition of a periodic potential. Numerically at least, finding the exact ground-state wavefunction 
is then a straight-forward application of one's favorite method — imaginary time propagation, shooting or relaxation. 
Similarly, propagation problems may be solved by an appropriate split-step algorithm for example. However, these 
procedures are computationally very intensive as compared to the analogous problems in the absence of the lattice for 
the following reason. Typically, we are interested in cases for which the period or "lattice constant" is relatively short 
compared to the extent of the wave function when just the trapping potential is present. (If that is not the case, one 
is really not dealing with a lattice at all, but a small set of coupled non-identical wells.) Thus we are concerned with 
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the case of a rapidly varying periodic potential superposed on the slowly-varying background of the trap potential. 
Hence as the wave function samples many wells of the potential, it possesses a detailed structure on the same scale as 
the lattice and a large number of grid points is necessary for accurate calculations. A time-dependent propagation of 
the field in a two-dimensional potential say, represents a formidable calculation. While such calculations are certainly 
possible JH , it is desirable to find a method that would circumvent the computational burden as well as highlight the 
underlying physics of the system. 

Consider first the stationary problem, where it is clear what the general form of the ground-state solution must be: 
to minimize the energy, the atoms pile up in the lattice wells, and the relative density in each well is determined by 
the local value of the confining part of the potential. The natural way to treat such a problem is to capture the rapid 
variation on the scale of the lattice using the Bloch functions of the periodic potential. If the nonlinearity is not too 
strong, the solution may be represented as a sum of a small number of slowly-varying envelope functions modulating 
appropriate Bloch functions. These envelope functions obey simpler equations from which all the rapid variation has 
been removed. For propagating problems, one may proceed in the same fashion but using different Bloch functions. 
In the language of solid state physics, this approach is essentially an effective mass representation with the addition 
of the nonlinear perturbation of the atomic collisions. We should remark that the work of Zobay et al. j| also uses an 
envelope function approach, but their work does not include a confining potential and is restricted to solutions close 
to a band edge. In addition, they use plane waves rather than Bloch functions as a basis, which limits the treatment 
to only relatively weak lattices. 

The paper is structured as follows. In section Owe derive simple approximate equations for stationary and propa- 



gating BECs in a finite optical lattice. In section [II we present solutions for the stationary ground state, comparing 



exact numerical solutions with our Bloch function approach. Section IV discusses extensions to treat fields near a 
band gap before we conclude. 

II. BEC IN A FINITE OPTICAL LATTICE 

We consider a BEC in an infinite optical lattice with an additional confining potential. To avoid incoherent heating 
from spontaneous emission, the laser beams creating the lattice are detuned far from the atomic resonance. In this 
case, the atomic field for the internal excited state may be adiabatically eliminated and the optical lattice acts as 

an external potential for the ground state field V> ff (x) JnjjljJ . Then the macroscopic wavefunction ^>(x) = /^ 9 (x) 

satisfies the CPE 
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V 2 + E/(x) + V(x) + I>(x, t)f 

2m 



V>(x,t), (1) 



where the optical dipole potential U(x) is periodic in one or more dimensions. For simplicity we assume that the 
beams are aligned along the coordinate axes and have wave vectors k x , k y and k z . Thus dropping an unimportant 
constant in the potential we may write 

U(x) = -^KpCos(fc p iz;p) (2) 

where runs over the coordinates x,y, and z, and fc M = 2fc M = 27r/d jU , where is the period of the lattice along 
axis [i. Also, V^(x) is the magnetic trap potential and T = 4arK a/m is the interatomic interaction strength, with a 
the s-wave scattering length. 

We now seek solutions to Eq. (|l|) using an expansion in terms of suitable Bloch functions of the optical potential. 
It will be helpful to begin by stating a few results. 

A. Band structure and Bloch functions 

We consider stationary solutions to the linear Schrodinger equation 

h 2 



--^V 2 + C/(x) 
2m 



x)=^(x). (3) 
Bloch's theorem states that the solutions take the form 



x) = </Vk(x) = e lkx u„ k (x), (4) 
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where the functions it„k(x) = u„k(x + R s ) share the periodicity of the lattice potential U, ie., R s is any "lattice 
vector" , translation by which maps the potential onto itself. The energy eigenvalues E = E^k form a series of bands 
indexed by the quantum number n, with an infinite number of bands for each value of the wavevector k. It is also shown 
in any solid-state physics textbook, that the wavevectors are unique only up to the addition of an arbitrary "reciprocal 
lattice vector". One can then restrict consideration to wavevectors lying within the "first Brillioun zone" Jin . In 
the one-dimensional case, this corresponds to the range \k\ < ir/d. We shall use a Dirac notation in which the Bloch 
function (/>„k(x) is denoted by |nk). In the usual way we restrict the wavevector to N discrete values with N large, 
so that the Bloch functions are normalized as 

(mk\nk'} = / rf 3 x^ k (x)^ k /(x) = ATW mrl c>kk' , (5) 

JN colls 

where the integral is taken over TV unit cells of volume f2. 

For the potential in Eq. @, Eq. (|) is separable and the well-known solutions are given by Mathieu functions [Q_3| . 
In practice, the solutions and eigenvalues are probably found most simply by solving Eq. (ra) directly. 



B. Equations of motion 

1. Bloch function expansion 

Using the properties of the Bloch functions, we now develop an envelope function approach for an atomic field 
ip(n, t = 0) in a finite lattice. Certainly any wave function ?/>(x, £) may be expanded in terms of the complete Bloch 
modes as 

^(x, *)=££ A nk (t) |nk) e - iE ^' n . (6) 

n k 

One could easily write down evolution equations for the amplitudes A nk but such a procedure is of little use for finite 
systems in which the wave function does not have the infinite extent of the Bloch modes. A very large number of the 
A n k would be required to accurately describe the dynamics. Instead it is more natural to use the Bloch functions to 
capture the rapid oscillations in the wave function but describe the localization of the field through slowly varying 
envelopes. We thus suppose that the atomic field is characterized by a central wave vector ko and try an expansion 
of the form 

V(x) = ^/ nko (x,i)|nko)e- £ -o*M. (7) 

n 

where the envelope functions /„k (x,i) are assumed to vary slowly on the scale of the lattice period. This is just the 
case for atomic BECs or an atom laser, as the momentum distribution is always distributed narrowly about some 
central wave vector. We defer the question of the choice of ko and the calculation of the initial envelopes / n k (x, t = 0) 
for a moment and concentrate on the dynamics. 

Substituting Eq. (fy into the governing Eq. |l]) gives a complicated equation involving sums over all the envelope 
functions and still involving the rapidly oscillating periodic potential. To obtain a simpler description, one insists 
that the variation of the /„k with x occur on slower scales than that of the Bloch functions |nko). Such a condition 
may be applied using the effective mass formalism of solid state theory |l2|, or equivalently, using "multiple scales" 
techniques as has been demonstrated in the context of nonlinear optics in a series of papers by de Sterke, Sipe and 
co-workers [p|Jl^|. Essentially, one obtains an equation for each envelope by dotting from the left with each Bloch 
function (nko|. 

In fact, if more than two envelopes are included the resulting equations remain quite complicated. For a start, all 
the bands are coupled by the nonlinearity, but there are also linear couplings that occur between bands with a small 
energy separation p3|. However, for many problems of interest, it may be sufficient to consider just one or two bands. 
For example, the stationary ground state solution must attain the lowest energy possible. Providing the nonlinearity 
is not too large, the dominant contribution to the solution is then provided by the lowest band at ko = 0. If the 
nonlinearity is very large, coupling to the higher bands can not be ignored and either more bands must be included or 



a different approach adopted. We return to this point in section IIIB. A one band approach can also be acceptable 
for time-dependent problems. A natural experiment might be to form a condensate in the stationary ground state of 
the lattice and then apply a weak momentum kick via a momentary change in the potential. In this case, we would 
expect the solution to be well represented by an expansion on the lowest band but concentrated around a nonzero 
value ko. If the system is tuned close to a band gap, however, there is strong coupling between the Bloch modes 



above and below the gap, and a two- mode approach is necessary M. We consider this situation in section IV 
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2. Single band approach 



In calculating solutions we shall concentrate on the stationary ground state of the BEC in the finite lattice, and 
hence a single envelope function is sufficient. We emphasize that for the envelope function approach to be valid, the 
nonlinearity must be considered a pertu rbation. We discuss the regime of applicability imposed by this constraint 
when we examine results in section III A . 

By applying the effective mass or multiple scales methods we simply obtained the reduced GPE for the 

envelope function / n k j 



iti 



d/nkp 

dt 



V/„ ko = - 



if 



2m*^ dxndxu 



/nk + ^( x )/nk + T* |/„ ko | 2 /nk , 



where 



v 9 = v nn (ko), 



f 



1 | ^(ko)< p (ko) + ^„(k )< p (ko) 

pjtn 



(E n k - Epk ) 

r* = £/ d 3 x| Unko (x)| 4 , 

" J cell 

and the velocity matrix elements are defined as 



v y( k o) 



(ik |p| jk ) _ l_ 
m 



m 



TibaSij + / d 3 xM iko (x)*p!i jko (x) 

./cell 



(8) 

(9a) 
(9b) 

(9c) 
(10) 



Equation (|^) is valid provided that the maximum nonlinear "detuning" r*|/„k | 2 a t the point where |/ n k | attains 
its largest value is small compared to the energy separation to the nearest adjacent band at wave vector ko. If this is 
not the case we must use a two band model |qji"4[ as discussed in section IV. We now have a significant advance in 
that the periodic potential does not appear explicitly in Eq. (|^) . Its action is effectively represented through the three 
parameters in Eqs. (^) which we now discuss. The envelope /„k moves with a group velocity v s which is simply the 
gradient of the band n at wavevector k . Note that at particular points of symmetry the group velocity vanishes, 
notably at the band edges, and of course at the bottom of the lowest band for n = 0, k = 0. The bare atomic mass 
m has been replaced by the effective mass tensor m* which is determined by the interaction of the other bands with 
the band of interest [Eq. (pq)]. The term involving the effective mass in Eq. (||) acts as a source of dispersion on the 
atomic field. Taken together the group velocity and effective mass constitute a quadratic expansion of the band n 
around the point ko. Finally, the parameter T* describes the change in the effective nonlinearity due to the shape 
of the Bloch function. With a weak potential, for which u„k is almost a constant function, T* w T. However with 
a strong lattice, for which the Bloch functions become peaked within or between the wells, the magnitude of the 
effective nonlinearity increases significantly. 

The wavevector ko is chosen in the first Brillioun zone and may be known from the nature of the problem — say an 
atom laser pulse of known momentum incident on the optical lattice — or might be calculated using an appropriate 
definition such as 



k 



dk 3 |i/i(k,i = 0)| 2 K(k) 



(11) 



where ipik) is the Fourier transform of '(/'(x) and K(k) is the vector in the first Brillioun zone equivalent to k. 
stationary state is sought we of course have k = 0. 



If a 



3. Soliton solutions 

Equation (||) is a nonlinear Schrodinger equation (NLSE) describing the evolution of the envelope field as it prop- 
agates through the lattice under the influence of the atomic interactions and the external potential V(x). Note that 
for one-dimensional geometries in the limit of vanishing V(x) we can obtain soliton solutions for the envelope. It 
must be emphasized that it is the envelope function that has soliton solutions, hence these are not the same objects 
as the standard atomic solitons found when the field operator or wave function itself satisfies a NLSE Jl5| — . The 
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envelope soliton propagates at a velocity determined by the local gradient of the band structure and with a width 
determined by a balance between the nonlinearity and the effective mass, rather than the bare mass. Such "Bragg 
grating solitons" have been observed in recent years in optical fiber gratings at extremely high laser powers [p|.^9[. 
In the present case, the nonlinearity is much stronger and BEC grating solitons should be easier to generate. For 
example, our calcul ations suggest that for sodium atoms in a weak lattice generated by a diode laser of wavelength 
985 nm (see section [II A ), a condensate of 500 atoms could produce a Bragg soliton extending over 30 periods of the 
lattice. 

For time-dependent problems, the term in V(x) in Eq. Q describes the response of the soliton to slow changes in 
the non-periodic potential. For example a broad laser beam incident from the side could be used to produce a local 
"hill" in the optical lattice. Gradients in this potential of course generate spatial oscillations in the field, which is to 
say, they induce wave vector shifts. For modest wave vector shifts, Eq. (^) correctly describes the resulting changes in 
velocity and shape of the soliton. However, if the shift becomes too large, the quadratic expansion of the Bloch band 
implicit in the effective mass approximation breaks down, and there must be a dynamic reallocation of the coefficients 
v g , m* and T* corresponding to the new point on the band |2(J. Alternatively, one could extend the effective mass 
approach by including higher spatial derivatives in Eq. (|J). 



III. STATIONARY SOLUTIONS 



While the possibility of Bragg solitons in BEC's is interesting, it is the properties of the ground state solution of a 
BEC in a lattice which are initially likely to be of most importance for experiments, and it is on these that we focus. 
We shall concentrate on the one-dimensional geometry which is perhaps the most pertinent experimentally. The 
harmonic trap potential I^(x) is assumed to be tightly confining in the directions transverse to the lattice potential 
and weakly confining along the lattice, such that the trap angular frequencies in the perpendicular (uj r ) and axial 
(u! z ) direction satisfy uj t ^> ui z . Eq. (0) then reduces to 



h 2 d 2 
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r^r - ncos{k z z) + -mcuiz^ + NT\ip(z)\ 

2m dz £ 2 



(12) 



where T = T/A e g = 2huj r a is the effective nonlinearity taking account of the transverse mode profile. There are two 
natural choices of dimensionless variables: the usual set of harmonic oscillator units and a set in which the period 
of the lattice is unity. We shall find it useful to use both at different times. The sets are i ntroduced by choosing 
appropriate energy and distance scales. For the oscillator case, these are e D = Tilo z and Z D — \Jti/ (muj z ) respectively. 
For the lattice case we instead have ei = h 2 /(md 2 ) and Z\ = d. We will use an overbar to denote quantities in 
oscillator units and a caret for quantities in lattice units. Thus in oscillator units we have £ = z/Z , /(C) = ^fZ^%l)(z), 
A = XZ , R, — n/e , £l — /i/e , = ?iw 2 /e = 1 and C — NT/(e Z ). Similarly for the lattice case, we have 
£ = z/Zi, g(£) = ^fZ~iip(z) and A = XZi = 2n, while the other parameters are defined in the obvious fashion, replacing 
the subscripts o by I and the overbars by carets. 

We have included the number of atoms N in the nonlinear constants C and C and so the wave functions are 



normalized to unity. We now turn to finding the ground state solutions of Eq. (12), identifying a number of regimes 
depending on the relative strength of the interactions and the lattice. 



A. Weak nonlinearity 

The regime of a weak nonlinearity is handled well by the effective mass envelope function approach introduced in 
section pi In this case the nonlinearity and the harmonic trap both act as perturbations to the dominant lattice 
potential and it is convenient to use lattice units. Thus we seek solutions to 
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l d 
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_fccos(27rO + 5 n^ + C| S (0| 



<?(0- 



(13) 



The wave function is written g(£) = a(£)0o(£) where 0o(£) is the lowest energy Bloch function (k = in the lowest 
band) for a lattice of depth k and period 1, and has eigenvalue Eoo(k) = £ , n= o,fc=o/e;- In the lattice units, Eq. (|^) 
becomes 



Ea(0 



1 



2m* d£ 2 



C*|a(0P 



a(0. 



(14) 
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As the lattice period in Eq. (|12|) is unity, the shape of the Bloch function </>o(£) is specified by the single parameter 
as are the values of the nonlinearity C* = j(k)C, with 



1 
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7 (fc)= dtlMOF, (15) 







and the (now scalar) effective mass given by 

1 =1 , 2 J2 l{MOHd/d ^ MO)l \ (16) 
m * p >o ^oo - E p0 

where the Bloch functions 4> P {£,) and energies E p q are of course normalized in lattice units. If we did not use lattice 
units, we would be forced to consider the dependence of these quantities on the lattice period as well. The effective 
quantities for the cosine potential are plotted in Fig. |asa function of k. Note that in the limit of a vanishing lattice 
potential, rh* — 1 and j(k) = 1. For large k, the Bloch function tends to a sequence of harmonic oscillator ground 
state wave functions and — * k x l A . 

The importance of Eq. (|14[) is that the problem including the lattice has now been reduced to the usual GPE for a 
condensate in just a harmonic trap with no lattice. All the physics involving the lattice is contained in the effective 
parameters and we are left with an equation whose properties are well understood. Equation ( p^ ) may naturally 
be solved in identical fashion to the standard GPE — numerically in general, or using Gaussian or Thomas-Fermi 
approximations (TFA) in appropriate limits. In particular, the TEA solution has the form 

K61 2 = max P-gf /2 ,o|. (17) 

Note that as the effective mass grows rapidly with k (see Fig [|), the TFA to Eq. (|IJ) is valid for quite modest 
values of k. It must be emphasized though, that physically, this "effective mass TFA" (EMTFA) is a rather different 
approximation to the standard TFA which applies when the total kinetic energy is negligible. Here it is only the kinetic 
energy associated with the envelope junction which may be considered small. Indeed there is considerable kinetic 
energy associated with the oscillation in the Bloch function but this of course has been removed from Eq. (|l4|) . It must 
also be realized that even if Eq. (0) is solved numerically, this represents a significant gain over a direct numerical 
solution of the complete GPE with the lattice potential included. From a practical point of view, the numerical routine 
is much faster as it requires a small fraction of the grid points required to simulate the whole lattice. This advantage 
is most obvious when one considers time-dependent problems using Eq. (^) rather than Eq. (|l|). Less prosaicly, the 
separation of scales allows a clearer physical insight into the expected general shape of the ground state condensate, 
with the precise details of the rapid oscillation subsumed into a few parameters, that need only be calculated once. 
Let us now consider the regime of validity for the effective mass approach as a whole and the EMTFA. To justify 

13| ) must be small 
r the nonlinearity 



an expansion in terms of the linear Bloch modes of the lattice, the effect of the nonlinearity in Eq. ( 
compared to that of the lattice potential, on the scale of a single period. This is the sense in whie 
is to be regarded as a perturbation. Thus if p = ,g(£) 2 is the density of the envelope function, we require 

dp < k. (is) 

Assuming that this is so, it may be checked that the EMTFA to Eq. ([l4|) is valid if 

C* » 2 V ^/(m*) 3 / 4 , (19) 
Using the EMTFA expression for the peak envelope density at the origin, Eq. (O) may then be expressed as 



^oC« -V7(«0(2£) 3 . (20) 

Note that due to its weak dependence on k, the dependence on ^(k) may usually be ignored. If the EMTFA is not 
valid, we can obtain a different bound on C by noting that the peak density of the envelope will certainly be lower 



than for the associated linear system for which a(£) = 1/ ^to^/tt exp(— £ 2 / (2rg)) , with r = l/ytlom*. Thus we 
obtain the condition 



£<W— • (21) 
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But as Eq. ( |l9| ) is not true, Eq. (|T]) is satisfied provided only that *C \pHk. These conditions are summarized 
in Table | along with the chemical potential which is given by /i = E 00 {k) + E -> £ 00 + (3C , *J7 /2) 2/3 /2, where the 
second relation holds when the EMTFA is valid. To demonstrate that these conditions could be satisfied by current 
experiments we consider Na 23 atoms with a lattice produced by counter propagating beams from a diode laser at 
985 nm pl| so that d = 0.49 /im. Taking a detuning of A = 10 10 s~l and a maximum Rabi frequency of O e = 0.01 A 
we obtain a maximum lattice depth of k max = K max /ei = fr^e/( e iA) » 90. Assuming w r /(27r.50) = lo z /(2tt) = 20 Hz 
and a scattering length a = 5 nm gives fio = 0-01 and C = 0.0137V. Choosing k = 90, we could satisfy Eq. ( po"| ) with 
condensates of up to a few million atoms. 

We show two examples of the envelope approach in Fig. |], comparing the condensate density as a function of 
position according to the various models discussed above. For both plots the nonlinearity C = 1 and trap frequency 
f?o = 0.05. Fig. ||a shows a case for a weak lattice with k = 10, while in Fig. ^b, the lattice strength is k — 100. 
The results of exact numerical calculations using Eq. (|l3|) are shown by the solid line in each figure. For clarity, we 
show only the positive £ part of the solution; the negative part is an exact mirror of this. The overall extent of the 
wave function is similar in the weak and strong lattices with the weak lattice case appearing to have a longer tail in 
the condensate density. In the stronger lattice, the oscillations are clearly more pronounced with the density confined 
more closely to the center of each well. 

Turning to the accuracy of our approximate descriptions, for the weak case in Fig. |^a, the effective mass description 
of Eq. (|l4|) is shown with a dotted line which is barely distinguishable from the exact result. The dashed line indicates 
the EMTFA solution which provides reasonable agreement with the exact solution near the center of the trap but 
fails towards the edge where the neglect of the envelope's kinetic energy with respect to the nonlinear part is invalid. 
For the strong lattice (Fig. the effective mass and EMTFA results are now perfectly superposed and are both 
indicated by the dashed line which has excellent agreement with the exact case for all but the last two peaks of 
the distribution. These results concur with the validity regime calculations above: due to the difference in effective 
masses, Eq. ([l9|) predicts that the EMTFA should be valid for the strong lattice but not for the weak. 

An additional dot-dashed line in each figure indicates the exact solution of Eq. ([l2]) when the nonlinearity vanishes 
((7 = 0). Comparing the two linear solutions shows how in the absence of a nonlinearity, the condensate shrinks with 
increasing lattice strength due to the increase in effective mass. In both figures, the exact solution is clearly broader 
than the linear result, despite the repulsive effect of the kinetic energy being reduced through the effective mass. This 
reduction is of course outweighed by the dominant repulsive nature of the nonlinearity. In spite of this, the linear 
Bloch function is perfectly adequate to represent the rapid oscillation — the nonlinearity is not strong enough to distort 
the wave function on the scale of a single period. 

Finally we consider the relative density distributions for the two exact solutions. Increasing the lattice strength 
produces two competing effects — the increase effective mass reduces the repulsion due to the kinetic energy, while the 
increase in effective nonlinearity increases the interatomic repulsion. For large enough effective mass, however, the 
kinetic energy plays essentially no role and any further increase in the mass makes no difference, whereas the effective 
nonlinearity and therefore the nonlincarly-induced repulsion grows slowly but without bound. Hence, we should 
expect that the condensate spreads out with increasing lattice strength. In Fig. [| we plot the cumulative density 
J_£ d£' \g(C)\ 2 as a function of £ for several values of the lattice strength k. From this plot it is clear that in spite 
of the weakening kinetic energy, the condensate slowly expands as the lattice strength ranges over k = (dot-dashed 
line), 10 (solid), 50 (dotted) and 100 (dashed), due to the increase in the effective nonlinearity [see Eq. ((L5|)]. The 
second dot-dashed line indicates the cumulative density according to the EMTFA envelope for k = 100. Note that 
the exact solution in the weak lattice appears to have a longer tail than for the strong lattice (solid lines in Figs. ||). 
Figure || shows that the structure of the tail is not indicative of the global density distribution. 

B. Strong nonlinearity 

As the nonlinearity increases or the trap tightens, the assumption that the nonlinearity plays no role over a single 
period [Eq. Jig)] breaks down. While one may use additional envelope functions each with their own Schrodinger 
equation, the utility of such an approach decreases quickly. However, if the nonlinearity is sufficiently large, we may 
abandon the Bloch function approach altogether, instead using the standard Thomas-Fermi approximation (TEA) by 
neglecting the kinetic energy term in Eq. ( |l2"| ) . S0renberg and M0lmer have used this approximation for the infinite 
uniform lattice. Here, we examine the more complicated system including the trap potential in the Thomas-Fermi 
regime. In this regime, the lattice is effectively a perturbation to the dominant harmonic trap, and it is natural to 
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invoke oscillator units, so that Eq. ( |12| ) becomes 

Id 2 I 



M/(C) 



.___ Scos( ac) + -c 2 + c|/(C)| 2 



/(C)- (22) 



|/(C)| = max | 2 = ,0|. (23) 

In estimating the boundary point Co of the condensate we may neglect the effect of the lattice, and take the usual value 
Co = y/^fi- This is consistent with the nature of the approximation which is poor at the boundary anyway. Moreover, 
to the same degree of approximation the dependence of the chemical potential on atom number is unchanged from 
the lattice free result fi = (3C/2) 2 / 3 /2. 

The regime of applicability for the TFA is analyzed as follows. Wc must certainly have fi > R which prevents the 
solution from attaining zero density regions in every lattice well. The standard TFA condition that the total kinetic 
energy be negligible compared to the total potential and nonlinear energies is given by (7 > 1 as for the lattice- 
free case. However, if the oscillatory part of the solution is to be accurately represented by the factor Rcos(A()/C 
in Eq. (^3|) we must also consider the relative energies associated with just the oscillatory part. An approximate 
condition may be found by evaluating these energies over a single lattice period at the center of the trap and insisting 
as usual that the resulting kinetic energy be small compared to the potential and nonlinear energies. Assuming R -C fi 
we eventually arrive at the additional condition 

fi » A 2 /4, (24) 

or 

C > 0.24A 3 . (25) 

Note that in the lattice units, this has the slightly odd form fi ^> it 2 . We see below that this new condition is genuinely 
necessary. 

We show examples in Fig. [| with exact numerical solutions to Eq. (|2^) in solid lines, and the TFA in dotted lines. 
In all three plots we have C = 3000. In Fig. ||a the lattice strength R = 10, and A = 2ir. The exact chemical potential 
is fi = 136.2. The requirements for the TFA are thus satisfied and there is good agreement between the exact and 
approximate solutions except near the boundary of the condensate. Figure |i|b shows the solutions where the lattice 
strength is increased to R = 50, for which the exact chemical potential is fi = 133.8. The region of disagreement 
near the boundary is now of course larger, but there is still good agreement near the center of the condensate despite 
R approaching 50 % of the chemical potential. In Fig. [|c, we have R = 50 with a reduced lattice period such that 
A = 57T. The chemical potential is fi = 134. As the ratio fi/(A 2 /4) = 0.54, Eq. (Q) indicates the TFA is invalid, 
in accord with the poor agreement between the TFA and exact curves. We emphasize that if the lattice were not 
present, the TFA would certainly be valid in this case because C » 1. Hence the condition given by Eq. (|24|) is 
indeed important. The dashed lines in Fig. ||c indicate the extent of oscillations in the solution predicted by the single 
band Bloch function description used in the previous section for weak nonlinearity. This is clearly unsuccessful in this 
regime as we would expect. Comparing it to the exact solution indicates that the repulsive nature of the nonlinearity 
tends to reduce the amplitude of the oscillations as compared to the linear regime. Note that despite the range in 
shapes of the three exact solutions in Figs, ^a-c, the chemical potential varies by only a few per cent from the lattice 
free value fi = 136.2. The parameter values used here would be easily obtained experimentally. Using the same trap 
frequencies w r /(27r.50) = lo z /(2tt) = 20 Hz, we find a maximum lattice strength of R m ax ~ 8000 and an effective 
nonlinearity of C = 0.1N. Hence, a condensate of about N — 30000 atoms would correspond to the predicti ons in 
Figs. [|. Note that when expressed in equivalent units, the lattice strength here is around 1 % of that in section III A, 
which explains why condensates of a million atoms could be considered to be weakly nonlinear in that section. 

It is worth noting that in the TFA regime, an envelope function approach is not appropriate, because the amplitude 
of the oscillations in the condensate density due to the lattice are of the same size across the whole condensate. In 
contrast, we expect an envelope function approach to be appropriate if the rapid oscillations in the density scale 
linearly with the envelope density. This difference is clearly apparent if the exact solutions in Figs. || and ^ are 
compared. 
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IV. DISCUSSION: ATOMIC FIELDS NEAR A BAND EDGE 



Finally, we return to the propagation problem and briefly discuss the rather different physics involved for condensates 
with a center wave vector tuned close to a band gap. In this case, unless the gap is very large and the field is tuned 
close to one edge, then neither of the two bands surrounding the gap can be regarded as remote and the field should 
be expanded using a Bloch function from each band. As our primary interest has been in the ground state stationary 
solution which is always far from any gap, we shall not explore this case in detail. For the sake of completeness, 
however, we state the dynamical equations that correctly describe such a situation using a method of de Sterke et 
al. @. 

We label the two bands u (upper) and I (lower) and expand around a frequency u>$ = (uj u k + cj; k )/2: 

^(x, t) = [/ uk (x, t) \vk) + /, k (x, t) |Zk)] e-^K (26) 

Again one insists that both / uk and /; k are slowly-varying functions and using the effective mass or multiple scales 
techniques obtains two coupled nonlinear equations for / uk and // k . In coupled mode theory of periodic structures it 
is often found convenient to work with envelopes modulating forward and backward going amplitudes as opposed to 
Bloch functions as these envelopes correspond closely to the field incoming and outgoing from the periodic structure. 
It can be shown that the amplitudes G± = (/; k =F i/uk)/2 correspond to such amplitudes. In the limit of a weak 
potential where the Bloch functions at the band edge are cosine functions, the amplitudes G± simply modulate forward 
and backward-going plane waves. In terms of these variables the final dynamical equations are [|l4| 



dG 

i 



dt 
dG. 



v 3 • VG+ + <rG_ + V{x)G+ + r (|G+| 2 + 2|G_| 2 )G+ 
- v g • VG_ + ctG+ + y(x)G_ + r (|G_ | 2 + 2|G+| 2 )G 



FiGG+l 2 + \G- | 2 )G_ + Y X {G+G*_ + G;G_)G + + T 2 G 2 _G* + = 0, (27a) 



-l-rxOG+l 2 + |G_| 2 )G+ + T^G+Gl + G;G_)G_ + T 2 G\G*_ = 0. (27b) 

The linear coupling parameter a — A/2, the group velocity v g = \ (uk|p| Zk) |, while the nonlinear coefficients are 
found from overlaps of the (real) Bloch functions: 

1 = ^ , (28) 

Qllll ~~ Q-uuuu / r)n \ 

i l = 2 ' { ' 

^ ocuii — 6a uu u + a uuuu 

1 2 - ^ , (oUj 



where 



/ d 3 X(/>i(x)0j(x 

icell 



Oijkl = q / ' / ' i x";('xK> ( -lx').v,lx)o/|x) (•'!!) 



for {i,j, k, 1} in {u, I}. In these equations, we have dropped the kinetic energy terms involving the effective mass 
altogether. The coupling between the bands through the parameter a itself produces a dispersive behavior which is 
much larger than the intrinsic dispersion within each band. In the limit of a weak potential, the coefficients Ti and 
T 2 vanish and one obtains the well-known coupled mode equations for propagation in nonlinear periodic systems. 
Such equations were recently obtained in the BEC context by Zobay et al. || and used to predict the occurrence of 
gap solitons in condensates with multiple internal levels. Their system differs from ours in that it involves a Raman 
coupling between two magnetic ground states whereas we have used a strictly scalar field interacting with an external 
potential. The fact that both systems can be reduced to identical equations is an indication of the underlying physical 
consequences of combining a nonlinearity with a periodic system — properties such as gap solitons and switching 
appear ubiquitously. The full "deep grating" equations J27| ) with r!,r 2 ^ (but V^(x) = were first derived by 
Salinas et al. |^2|| to describe optical propagation in Bragg gratings with large index variation between the two media. 
The additional terms in these equations lead to predictions that depart from the shallow grating results for index 
variations of 10 % or more. In the atomic case, the analogous situation is when the lattice potential is deep and 
the Bloch functions at the band edge depart from simple plane waves. The exposure of condensates to such deep 
lattices is certainly a realistic goal in the laboratory and the use of the deep grating equations may thus be important 
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in modeling propagation in lattices in the future. These arguments hold equally for the Raman geometry used by 
Zobay et al. 

We should remark that the use of a two band model as described here is not always possible. In the case of extremely 
deep potentials the band gap widths may broaden so much they become comparable to the original separations between 
the gaps. In that case, the existence of two bands that are significantly closer to the center frequency than all the 
others is doubtful. However, if the desired center frequency is closer to one band edge than the other, a one band 
model may again become feasible. 



V. CONCLUSION 



We have introduced a Bloch function analysis of a BEC trapped in a finite optical potential. The approach describes 
both pulse propagation and stationary solutions providing a single Schrodinger like equation for the envelope function 
provided the nonlinearity is not too large. If the physical parameters lie within the constraints we have given, the 
envelope function solution agrees extremely well with the exact numerical results. For very strong nonlinearities, the 
standard Thomas- Fermi approximation is the correct description for the ground state, but its application requires a 
stronger condition than in the absence of the lattice. 
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FIG. 1. Effective nonlinearity V (solid) and mass m* (dotted) for a one-dimensional cosine potential of strength k. 
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FIG. 2. Wave function as a function of £ for a) k = 10 and b) ft = 100. Other parameters are flo = 0.05 and (7 = 1. 
Line styles indicate models as exact numerical (solid), effective mass (dotted), EMTFA (dashed) and linear (dot-dashed). In 
b), the effective mass and EMTFA models are both represented by the dashed curve. 

FIG. 3. Cumulative density d£' \g(£,')\ 2 as a function of £ for values of the lattice depth k = 10 (solid line), 50 (dotted) 
and 100 (dashed). The dot-dashed lines indicate the cumulative density of the exact solution for k = 0, and of the EMTFA 
envelope for k = 100. 

FIG. 4. Exact wave function f(() (solid line) and Thomas-Fermi estimate (dotted) with nonlinearity C = 3000 and a) 
k = 10, A = 2tt, b) R = 50, A = 2-k and c) R = 50, A = 5-7T . The dashed line in c) denotes the prediction of the envelope 
function approach. 



Regime 


Conditions 


Chemical potential 


Weak: Cp < k , (EMTFA) 


C* » 2^n /(m*Y'\ QoC « | V7(^)(2£) 3 


p = £ 00 (k) + (3C" /2) 2 / 3 /2 


Weak: Cp < k , (non-EMTFA) 


C < kyj&orh* /tt 


A = E o{R)+E 


Thomas- Fermi : 


p > A 2 /4, C > 0.24A 3 


p w (3C/2) 2 / 3 /2 



TABLE I. Validity criteria and chemical potential for ground state solutions in different regimes. 
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Steel and Zhang, Figure 3 




Steel and Zhang, Figure 1 




Steel and Zhang, Figure 4a 




Steel and Zhang, Figure 4c 




Steel and Zhang, Figure 2a 




Steel and Zhang, Figure 2b 




